% Zernike多项式计算函数，输入最大项数j_max,圆孔直径D，分辨率N
% 输出ZM为泽尼克多项式矩阵，其第三维对应不同阶数的znk多项式
function ZM = ZernikeCal(j_max, D, N)

R = D/2;
x = linspace(-D/2,D/2,N);
[X,Y] = meshgrid(x,x);
[theta,r] = cart2pol(X,Y);       %转化为极坐标系
mask = (r<=R);        %限定计算范围是半径为D/2的圆
theta(~mask) = 0; r(~mask) = 0;
r = r/R; % 归一化至单位圆内
ZM = angular().*Barmak(); % Output Zernike Polynomials Matrix

    % 内嵌子函数,用迭代方法计算泽尼克径向多项式
    % % References: Recursive formula to compute Zernike radial polynomials, OL
    % Barmak Honarvar and Raveendran Paramesran
    function RM = Barmak()
        %RM(1,1) = 1;
        RM = zeros(size(r,1),size(r,2),j_max);
        RM(:,:,1) = double(mask); m_pre = 0;
        
        if j_max > 1
            for j = 2:j_max
                [n,m] = Noll_j_to_nm(j); %唯一对应,但从[n,m]到j不唯一
                if m_pre == m % 如果当前m等于上一次角向索引m，则径向多项式相等
                    RM(:,:,j) = RM(:,:,j-1);
                else % 根据Barmak's paper，径向多项式的迭代计算分为三种情况
                    if m == 0
                        j1 = Noll_nm_to_j(n-2,m);
                        j2 = Noll_nm_to_j(n-1,m+1);
                        RM(:,:,j) = 2*r.*RM(:,:,j2) - RM(:,:,j1);
                    elseif m == n
                        j1 = Noll_nm_to_j(n-1,m-1);
                        RM(:,:,j) = r.*RM(:,:,j1);
                    else % m ~= n
                        j1 = Noll_nm_to_j(n-1,m-1);
                        j2 = Noll_nm_to_j(n-2,m);
                        j3 = Noll_nm_to_j(n-1,m+1);
                        RM(:,:,j) = r.*(RM(:,:,j1)+RM(:,:,j3))-RM(:,:,j2);
                    end
                end
                m_pre = m;
            % end of the for-loop
            end
        end
    % end of the sub-function
    end
        
    % 内嵌子函数,用于计算泽尼克角向多项式
    function TM = angular()
        TM = zeros(size(theta,1),size(theta,2),j_max);
        for j = 1:j_max
            %---计算j对应的阶数n和角向频率m---
            %---------------------------------
            [n, m] = Noll_j_to_nm(j);
            
            if m==0
                TM(:,:,j)=sqrt(n+1).*mask;
            end
            if m~=0 && rem(j,2)==0
                TM(:,:,j)=sqrt(2*n+2)*cos(theta*m);
            end
            if m~=0 && rem(j,2)==1
                TM(:,:,j)=sqrt(2*n+2)*sin(theta*m);
            end
        end
        
    % end of sub-function
    end
    
% end of function
end